Phase field approach for modeling intracellular dynamics 
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We introduce a phase field approach for diffusion inside and outside a closed cell with damping and 
with source terms at the interface. The method is compared to exact solutions (where possible) and 
the more traditional finite element method. It is shown to be very accurate, easy to implement and 
computationally inexpensive. We apply our method to a recently introduced model for chemotaxis 
by Rappel et al. [Biophys. J. 83, 1361 (2002)]. 
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PACS numbers: 02.60.Lj,02.70.-c,82.39.-k,87.17. Jj 

In dealing with free boundary problems, the so called 
phase field approach 2] appears as a method of choice. It 
has successfully been applied to various problems ranging 
from dendritic solidification Q , viscous fingering and 
crack propagation In the spirit of time-dependent 
Ginzburg-Landau models, the method avoids the track- 
ing of the interface by introducing an auxiliary field that 
locates the interface and whose dynamics is coupled to 
the other physical fields through an appropriate set of 
partial differential equations. In comparison to the more 
traditional boundary integral methods, the method is 
much simpler to implement numerically. 

In this Brief Report we introduce a phase field model 
for intracellular dynamics i.e. diffusion inside and outside 
a stationary, closed domain with source terms at the in- 
terface. We apply the method to a recently introduced 
model for the response of a Dictyostelium amoeba fol- 
lowing stimulation with the chemoattractant cAMP pj. 
In , due to the need to use a finite element method the 
numerical implementation of the model was limited to 
two space dimensions and the cells were treated as disks. 
As we will see below, the phase field method is capa- 
ble of faithfully capturing no-flux boundary conditions. 
Thus, it becomes feasible to investigate more realistic cell 
shapes in three dimensions. 

Before we introduce our approach, we like to point out 
some possible extensions of our methodology. The phase 
field approach can be easily modified to include prob- 
lems where the domain boundary is not stationary. For 
example, force generation on cell membranes, leading to 
shape changes, can be incorporated within the phase field 
approach in a straightforward manner. This would re- 
quire adding an additional equation for the phase field, 
but does not require the explicit calculation of a bound- 
ary. We should stress that attempting to model this type 
of problem using conventional techniques, with explicit 
boundary tracking, become quite cumbersome. 

Let us first introduce the most salient ingredients of 
our method. Our purpose is to describe the situation 
where some chemoattractant diffuses only inside the cell, 
i.e. its concentration c obeys the diffusion equation inside 



a closed domain: 



dc 
di 
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and satisfies zero- flux boundary conditions: 
n ■ Vc = 

As a phase field we take a function that takes the value 
<pi n = 1 inside the cell and 4> ou t = outside the do- 
main and varies smoothly across the interface. A possi- 
ble choice in one dimension for a cell between x = — a 
and x = a is: 
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tanh((a - \x\)/£) 



(1) 



The variable £ denotes the interface width. We then de- 
fine a field u that is to obey the equation: 



, du 



D(j)Vu 



(2) 



Our claim is now that inside the domain the field u 
behaves very similarly to c. Let us therefore first show, 
for simplicity in one dimension, that in the thin inter- 
face limit one recovers the no-flux boundary condition. 
Integration of (J2J over the interface yields: 



dt 
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-£) » 0. From this we deduce: 



i 



and thus, in the sharp interface limit £ — > 0, the reflec- 
tive boundary conditions are recovered. It is also clear 
that inside the domain where <j> is constant u satisfies the 
diffusion equation. 

Another interesting property of Eq. @ is that the 
quantity J (j>u dx, which can be interpreted as the to- 
tal concentration inside the cell, is conserved under the 
dynamics. Thus, even if the field u may become non- 
zero outside the cell (and indeed it does), the total con- 
centration inside the cell remains constant. In fact, the 
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asymptotic solution of Eq. Q (together with phase field 
Eq. QJ) is u{x) — A where the constant A is determined 
through 2a A = J dx(f>u(t = 0). The real asymptotic so- 
lution is of course u(x) — A — J°^ a dxu(t = 0)/(2a). 
Because 0=1 inside the domain and <f> = outside of it, 
the error £ in the asymptotic solution depends solely on 
the amount of concentration initially near the interfaces: 
£ w dx<jnt(t = 0)/2a. Since <f> is an antisymmet- 

ric function around the interface, the error is minimal if 
we take in the phase field approach an initial condition 
u(t = 0) which is locally symmetric around the boundary. 

This raises the following important point. While the 
solution u outside the cell is not of physical interest, it is 
essential to keep track of it. In practice we solve (0 where 
the phase field <fi exceeds a small threshold S (typically 
5~ 1CT 8 ). 

As a first test for our model we now solve Eq. @ nu- 
merically in two space dimensions. For the phase field 
we take here: 

0(r) = i + itanh((r„-r)/0 (3) 

where r is the radius in polar coordinates. In the case 
of a radially symmetric initial condition we can compare 
the field u to the analytic solution v which is expressed 
in terms of the Bessel functions 

v(r,t) = ^a„J„(A„r)e" A "* 

a 

where A„ is defined as the smallest number for which 
J n {X n a) — and a n is the projection of the initial con- 
dition on the set of Bessel functions. We also compare 
both fields u and v with the solution w of a finite ele- 
ment method obtained with MATLAB's PDE Toolbox 
(The Mathworks, Natick, MA). 

In Fig.^Ji we show the spatial profile of the phase-field 
together with the exact solution at a given time. It can 
be seen that the agreement is excellent, with a relative 
error of less than 1% (see insert). In Fig. we show 
the fields u,v and w at three different points. Again the 
agreement is excellent. 

In view of these promising results, let us now consider 
diffusion in presence of a production term and damp- 
ing. Including a source term at the interface is rela- 
tively easy. We add to the right hand side of © a term 
b(V(j)) 2 / J dx{V4>) 2 . The factor (V</>) 2 ensures that it 
only acts at the interface and the denominator is a nor- 
malization factor. To compensate for the production we 
also include a damping term of the form —fj,(f>u. The 
equation of motion then becomes: 

<j£± = V • \d^u\ - fjufm + b-^ft— (4) 
ot L J J dx(Vcj)) 2 

We have compared the phase field method with these 
two supplementary ingredients with the finite element 
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FIG. 1: Comparison of our phase field model with both ana- 
lytic solution and finite element method. The equations are 
solved on a grid of size 201 x 201 with Euler's method, with 
D = 250, r = 5, Ax = 0.1, and At = 5 ■ 10~ 6 . The ini- 
tial condition is c = 1 + cos(7rr/ro). (a) Spatial profile at 
respectively at t = 0.004, t = 0.01 and t = 0.02 of the phase 
field and analytic solution. Since the curves are hardly distin- 
guishable, we mark the analytic solution by small circles and 
plot the difference of the profiles at t = 0.01 in the inset. The 
relative error of the phase field is seen to be smaller than 1%, 
being maximal at the boundary, (b) Time evolution of the 
concentration at r = 5, r — 3.37 and r — 1.66. The analytic 
solution is again marked by small circles. 



method, again in the case of a two-dimensional circular 
cell with ro = 5. As can be seen in Fig. [2 the result is 
excellent. 

We now use our phase field method to solve a biolog- 
ical model for the response of a Dictyostelium amoeba 
following stimulation with the chemoattractant cAMP 
|1| . In recent experiments the establishment of an asym- 
metry within a few seconds after a rise of extracellular 
cAMP was demonstrated. The cAMP however diffuses 
rapidly around the cell and the applied signal is several 
orders of magnitude larger than the value required to 
elicit a response. This strongly suggests the presence of 
an inhibitory mechanism that suppresses responses at the 
back. 

In an abstract model for the initial response of the 
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FIG. 2: Comparison of our phase field model with finite ele- 
ment method with source and damping. We have taken the 
same system size, initial conditions and parameter values as 
in Fig. Q and now b — 1000 and /j, — 50. We show here the 
time evolution of the concentration at r — 5, r — 4.05 and 
r = 0.86. The curves are again hardly distinguishable. We 
plot the error u p f — Uf e for r = 5 (the worst point of Fig. la) 
as a function of time. As one can see in the insert, the error 
goes initially rapidly up to around 1.5 %, but decreases then 
to around 6 ■ 10~ 2 %. 

cell to the chemoattractant was proposed: it was sup- 
posed that the membrane can be characterized in terms 
of three states: quiescent (with density p q ), activated 
(with density p a ) and inhibited (with density p{). Ini- 
tially the entire membrane is in the quiescent state. As 
the cAMP reaches the front of the cell the membrane 
becomes activated at rate a[cAMP] and an inhibitor, in 
suggested to be cGMP, is produced (at rate a g p a ). 
The inhibitor diffuses toward the back of cell where it 
turns the membrane from quiescent to inhibited with 
rate /3 r [cGMP]. Furthermore both activated and inhib- 
ited state decay spontaneously to the quiescent state at 
rates 8 and j3t respectively. The equations for the mem- 
brane state variables are thus: 

d p 

= -acp q + /3fPi - (3 r gp q (5) 
— = ac Pq - Sp a (6) 

= ~PfPi+ PrgPq + Spa (7) 

The reactants cGMP and cAMP diffuse respectively 
inside and outside the cell. At the membrane they satisfy 
zero-flux boundary conditions. There is an source term 
for the cGMP field that accounts for the production of 
cGMP at the interface. Both cAMP and cGMP fields 
are damped at rates p c and p g . The phase- field method 
is thus well suited to solve the dynamical equations for 
the cAMP and cGMP concentrations. For the phase field 
corresponding to the cAMP (which diffuses outside the 
cell) we simply take the complement of <j) given by Eq. ||3J| : 
<p c = 1 — 0. The equations for the membrane variables are 



solved on all lattice sites where (V<^) 2 exceeds a certain 
threshold, namely 10~ 4 . 

We now present a comparison of the phase field ap- 
proach and the results obtained with a finite element 
method in A circular cell of diameter 10 pm is placed 
in a square domain of 30 pm by 30 pm. The diffusion 
constant of cAMP and cGMP was taken to be identical: 
D c = D g = 250/.tm 2 /s. The values of the other parame- 
ters can be read in the caption of FigOJ 

In order to mimic the asymmetric cAMP stimulus we 
maintain the cAMP concentration at a value well above 
threshold at the upper left corner of our domain. As 
expected from our earlier results the agreement of the 
cAMP fields, that solely diffuse around the cell, is excel- 
lent, see Fig. OK- We observe a slight discrepancy for the 
cGMP field (also Fig.|30, which grows with time. This 
is related to the larger difference (of around 10%) for the 
membrane variables, the production of cGMP being pro- 
portional to p a . The origin of the discrepancy might lie 
in the way the state variables are calculated: on a ring of 
finite width in the phase field model, and on 40 points on 
the interface in the finite element method. At any rate, 
since the model is of a quite abstract nature and since its 
predictions are only qualitative, a detailed investigation 
is beyond the scope of this paper. 

From a computational perspective, we note that in the 
phase field method the CPU time required is linear in 
the number of lattice sites whereas it is quadratic in the 
number of nodes in the finite element method. One thus 
expects the phase field method to be much faster. How- 
ever this naive result must be taken with some caution as 
in the finite element method the nodes are not distributed 
uniformly in space. To resolve some particular region in 
space with high accuracy we are thus obliged to take 
a small lattice spacing in the phase field method, such 
that the number of lattice points exceeds the number of 
nodes. Other parameters that will affect the accuracy 
of both methods are the time step At and the integra- 
tion algorithm. Again a detailed comparison is beyond 
the scope of the Report. In practice it turns out that 
the phase field method where the lattice spacing is equal 
to the minimal distance between two nodes of our finite 
element method is about one order of magnitude faster. 

Finally, we turn to the chcmotaxis model in three space 
dimensions. For the sake of simplicity we take a spherical 
cell, i.e. with a phase field like in Eq. © but where r 
is now the radius in spherical coordinates. We consider 
again a cell of radius ro = 5pm, in a box of dimensions 
30 x 30 x 30pm. Now the stimulus is applied by setting 
the cAMP initially well above threshold at one side of the 
box, here taken to be x = —15pm. As is illustrated in 
Fig. the phenomenology of the two dimensional case is 
reproduced here. As the cAMP front progresses toward 
the back of the cell, only the front of the cell is activated. 

In conclusion we have proposed a phase field model 
for intracellular dynamics. Our method is shown to 
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FIG. 3: Comparison of our phase field model with finite ele- 
ment method in the chemotaxis model, (a) cGMP concentra- 
tion at front (upper curves) and cAMP at back (lower curves) 
of the cell as a function of time (b) state variables p a at front 
and back of the cell as a function of time. In the finite element 
method the grid consists of 216 nodes inside the cell (of which 
40 on the interface) and 543 on the outside (of which again 
40 on the interface). The grid and also the mass and stiffness 
matrices were generated by MATLAB's PDE Toolbox (The 
Mathworks, Natick, MA) after which the equations of motion 
are integrated with a Fortran code. The time step is taken to 
be At — 2 ■ 10 -6 . In the phase field method the equations of 
motion are integrated on a 151 x 151 grid, the lattice spac- 
ing thus being Ax = 0.2. Here the time step is taken to be 
At = 10 -5 . We have taken the following parameter values: 
a = 4 , p f = 0.01, fa = 0.533, S = 0.1, ju c = 0, p g = 0.12, 
and a g = 60000. 




FIG. 4: The cAMP concentration, represented by gray levels 
(where white corresponds to large concentrations and black to 
low ones) is for visual simplicity shown at the back of the box 
and is similar to the cAMP field surrounding the cell. The 
activity density of the membrane is shown on the sphere, also 
represented by gray levels. It is seen that whereas the cAMP 
has reached the back of the cell, it has not become activated. 
The equations of motion are integrated on a 61 x 61 x 61 grid, 
the lattice spacing thus being Ax = 0.5. Here the time step 
is taken to be At — 10 -4 . These results have been obtained 
for the same parameters as in Fig. [3] Only the production 
term has been increased in order to account for the increase 
in membrane surface. 



variable, that evolves under the appropriate Ginzburg- 
Landau type of equation. We are currently working along 
these lines of research. 
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be very accurate, easy to implement and computation- 
ally inexpensive. Another advantage lies in its much 
greater flexibility with respect to other methods, like 
the finite element method. We revisited a chemotaxis 
model and, when considering three dimensional cells, we 
reach the same conclusions as in |l| for two-dimensional 
cells. Whereas we have restricted ourselves to circular 
and spherical domains, the extension to other geometri- 
cal forms poses no major problems, the only task being to 
generate a phase field <f>. Even better, our approach can 
easily be extended to deal with non-stationary bound- 
aries. This situation arises in a multitude of biological 
problems. For example, the Dictyostelium cells change 
their shape continuously during chemotaxis. Another 
example are shape transformations observed in vesicles 
||. In these cases, the phase field becomes a dynamic 
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